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ABSTRACT 

We present an algorithm to photometrically calibrate wide field optical imaging surveys, that simultaneously 
solves for the calibration parameters and relative stellar fluxes using overlapping observations. The algorithm 
decouples the problem of "relative" calibrations from that of "absolute" calibrations; the absolute calibration 
is reduced to determining a few numbers for the entire survey. We pay special attention to the spatial structure 
of the calibration errors, allowing one to isolate particular error modes in downstream analyses. Applying this 
to the Sloan Digital Sky Survey imaging data, we achieve ~ 1% relative calibration errors across 8500 deg 2 in 
griz; the errors are ~ 2% for the u band. These errors are dominated by unmodelled atmospheric variations 
at Apache Point Observatory. These calibrations, dubbed "ubercalibration", are now public with SDSS Data 
Release 6, and will be a part of subsequent SDSS data releases. 

Subject headings: methods: data analysis, techniques: photometric, catalogs, surveys 



1. INTRODUCTION 

A common challenge for all physics experiments is relat- 
ing a detector signal to the underlying physical quantity of 
interest. Astronomical imaging surveys are no exception; a 
CCD camera counts Analog to Digital units (ADU) in each 
pixel, a quantity that is (approximately) proportional to the 
number of incident photons. This relationship must be cal- 
ibrated to yield physical flux densities (erg cm _2 s _1 Hz _1 ). 
Key scientific programs of current and next generation imag- 
ing surveys demand ever more precise photometric calibra- 
tions. For example, wide field imaging surveys allow one 
to measure the clustering properties of galaxies (and there- 
fore, the underlying dark matter) on scales otherwise accessi- 
ble only in the Cosmic Microwave Background (CMB); com- 
paring the CMB at redshift z ~ 1000 with the relatively re- 
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cent Universe at z< 1 allows increasingly precise tests of our 
cosmological model. The first such measurements of clus- 
tering on gigaparsec scales and larger were re cently reported 
dPadmanabhan et all 120071: iBlake et al.ll2007l) . These results 
emphasize the need for accurate photometric calibrations over 
wide areas; the underlying clustering signal is a rapidly de- 
creasing function of scale, and could easily be overwhelmed 
by percent-level systematic errors in the photometric calibra- 
tion. A second example is reconstructing the structure of the 
Galaxy, using the photometric properties of different stellar 
populations. There have been a nu mber of efforts to do this 
with existing data (lJuric et al.l l2005). and it is a key scientific 
program for the next generation of imaging surveys. Finally, 
there is the general (and powerful) motivation that reducing 
systematic errors invariably reveals hitherto unseen details 
and avenues of enquiry. Leveraging the current and next gen- 
eration of imaging surveys to yield their maximum scientific 
potential requires revisiting the problem of photometric cal- 
ibration, moving beyo nd the simplifications currently made 
(Stubbs & Tonrv 2006). Several surveys are photometrically 
calibrated to a few percent; the challenge for the next gen- 
eration of surveys is to deliver < 1% calibrations over wide 
areas. 

Photometric calibration involves relating the output of a 
CCD to the physical flux received above the Earth's atmo- 
sphere. For wide-field imaging surveys, we separate this into 
two orthogonal problems - "relative" calibration, or the prob- 
lem of establishing a consistent photometric calibration (al- 
beit in possibly arbitrary flux units) across the entire survey 
region, and "absolute" calibration, which transforms the rela- 
tive calibrations into physical fluxes. This separation is useful 
since there exist a number of applications (such as the two 
discussed above) that are relatively tolerant of errors in the 
absolute calibration, but demand precise relative calibrations. 
Current calibration techniques, which usually involve com- 
paring observations to "calibrated standards", do not respect 
this distinction, making it difficult to control errors in the rel- 
ative calibration. Furthermore, calibrating off standard sys- 
tems normally involve relating different telescope and filter 
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systems, and are quickly limited by the accuracy with which 
these transformations can be measured. Accurate relative cal- 
ibrations would therefore only use data from the native tele- 
scope/ filter system, obviating the need for any such transfor- 
mations. 

A s econd separation, emphasized by Stub bs & Tonryl 
(2006) is to separate the "transfer function" of the telescope 
and detectors from that of the atmosphere. The telescope 
and detectors form an approximately closed system whose re- 
sponses can be (potentially) mapped out with exquisite preci- 
sion with laboratory equipment. The atmosphere, on the other 
hand, is an open, highly dynamical, system with a range of 
relevant time scales; the best one can do is to monitor it with 
limited precision. Although we agree with this separation in 
principle, applying it would go significantly beyond the scope 
of this paper. We therefore do not make this distinction in the 
analysis presented here, but we return to it at the end of this 
work. 

Techniques for relative photometric calibration have been 
applied to optical imaging in the past, although much more 
limited than the present work in the scope of eit her t he num- 
ber of objects or the field of view. Landoltl d 1983b and lLandoltl 
(1992) are widely recognized as describing one of the best- 
established "photometric systems". Landolt observed several 
hundred stars near the celestial equator with a photomultiplier 
tube on the Cerro Tololo 16-inch and the Cerro Tololo 1.5-m 
telescopes. Landolt achieved exceptional relative photometric 
calibrations in five broad optical bandpasses (Johnson-Kron- 
Cousins UBVRI). His data are accurate to 0.3% per obser- 
vation of each star, with an even better accuracy implied for 
those stars that have many observations. Unfortunately, the 
full benefit of the accuracy of this photometric system cannot 
be realized for other surveys due to the systematic uncertain- 
ties in transforming from Landolt's system responses to ob- 
servations on other telescopes using (typically) CCD photom- 
etry. There are some observations using exactly the Landolt 
system, most famously of supernova 1987 A, which made use 
of the otherwise-decomm issioned Cerro Tololo 16-inch tele- 
scope dBlanco et al.lll987l) . 

The other example of accurate relative optical photometry 
has come from the searches for massive compact halo ob- 
jects ( MACHOs) from microlensing events in dense star fields 
(e^g JUdalski et all Il992t lAlcock et alJl!993b lAubourg et alJ 
1993). MACHO events were detected from differencing im- 
ages taken with an identical instrument over timescales of sev- 
eral years. More recently, these same techniques have been 
used to detect the optical transits of planets. With a proper 
treatment o f the correlated n oise properties in the time series 
of images dPont et al.ll2006l) . it is possible to detect transits 
with peak depths of only ~1%. Note that the challenges here 
are different from the wide field imaging case considered in 
this work, since one is interested in differences in photometry 
of a single star. 

Cosmic Microwave Background (CMB) anisotropy exper- 
iments also demand very precise relative calibrations. This 
accuracy is obtained with repeat observations of the sky, 
and cross-linked scan patterns. The redundancy thus ob- 
tained allows one to simultaneously solve for the CMB tem- 
perature at a given direction on the sky and the detector 
calibration parameters. In this paper, we propose adopt- 
ing this technique as a new approach to calibrating optical 
imaging surveys - replacing the CMB temperature fluctua- 
tions with the magnitudes of stars. Note that as this in- 
volves comparing multiple observations, this is a differen- 



tial measurement, and therefore only yields a relative cali- 
bration. However, while the absolute calibration still must 
be obtained by comparison against standard stars, this is now 
applied uniformly across the entire survey region. These 
ideas are not new t o optical astronomy; precursors may b e 
found in th e work o f iMaddox et alJdl990l):lHonevcutll (jl992): 
FongetalJ (f 19921): iManfroidl d!993l); iFongetal. (1994); 



Glaze brook et al.l d 1994b . What is new to this work is both 



the (large angular) scales to which the method is applied and 

the accuracies obtained. 

The Sloan Digital Sky Survey dYork et al.ll2000l SDSS) is 
one of the most ambitious optical imaging and spectroscopic 
surveys undertaken to date. It has imaged a quarter of the 
sky in five optical bands, and has spectroscopically followed 
up more than a million of the detected objects. This makes 
the SDSS both a scientifically rich data set and an excellent 
proving ground for the next generation of surveys. Accord- 
ingly, our goal in this paper is to develop the idea above in 
the context of the photometric calibration of the Sloan Dig- 
ital Sky Survey (SDSS). We begin by recapitulating aspects 
of the SDSS essential to this algorithm in Sec. [2] Sec.[3]then 
presents the details of the algorithm. We then assess the per- 
formance of our calibrations with simulations of the SDSS; 
the results are in Sec. [4] We then present a recalibration of 
the entire SDSS imaging data in Sec. [5] Sec.|6]announces the 
release of this calibration to the public, and Sec.|7]concludes 
with a discussion of the features and limitations of this work, 
as well as its applicability to the next generation of imaging 
surveys. Although we focus on the SDSS, we phrase our dis- 
cussion in terms that allow adapting the methods described 
here to arbitrary imaging surveys. 

2. THE SDSS 



The Sloan Digital Sky Survey (lYork et alj20 00) is an ongo- 
ing effort to image approximately a quarter of the sky, and ob- 
tain spectra of approximately one million of the detected ob- 
jects. The imaging is car ried out by drif t-scanning the sky in 
photo metric conditions dlfogg et alfeOOll) . usi ng a 2.5m tele- 
scope (Gunn et al. 2006) in five bands (ugriz) (Fukugita et al. 
1996; ISmith et al.l 12001 using a specially designed wide- 
field camera dGunn et al.lll998l) . These data are processed 
by completely automated pipelines that detect and measure 
photometric pro perties of objects, and astrometrica lly cal- 
ibrate the data dLupton et aflbooU iPier et all 120031) . The 
first phase of the SD SS is complete and has produced seven 
major data releases (Stoughton et al. 2002; Abazaiia n et alJ 
[2001 [20041 l200l lAdelman-McCarthv et all I2006L 120071: 
Adelman-McC arfhv & for the SDSS Collaborationl2007l) lb . 

The SDSS imaging data (see also Fig.[T} are taken by drift- 
scanning along stripes centered on great circles on the sky in 
all five filters. These stripes are 2.5° wide, and are filled by 
two interleaved strips. The actual data is taken in runs, which 
are part of strips; multiple runs may be taken in a single night, 
not necessarily on the same strip. Each run is further subdi- 
vided into six camera columns or camcols, corresponding to 
the six columns of CCDs on the camera. The data from each 
CCD is in turn split into frames, consisting of 1361 drift scan 
rows. The five frames corresponding to the same region of 
sky observed in the five SDSS filters, are collectively referred 
to as afield. Note that while runs and camcols correspond 
to physical separations of the data, the division into frames is 
purely artificial. The integration time is approximately 54. 1 

16 http : / /www . sdss . org/dr 6 
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FIG. 1 . — The geometry of the SDSS imaging. Part of an SDSS stripe with 
the two interleaved strips (denoted by N and S) is shown. Each strip consists 
of six camcols (numbered 1 through 6 in the figure), while each camcol is 
further divided into fields (for simplicity, we show field divisions for only 
two camcols). See the text for more details. 



seconds per frame in each filter, with a time lag of ~ 73 sec- 
onds between each adjacent filter. The order of the filters as 
they observe the sky is riuzg. 

The current survey flux calibrations are applied in a three 
step process, involving three different telescopes and subtly 
different filter sys tems. The absolute flux system is defined 
by BD+17°4708 dOke & Gunnlll983h . an FO subdwarf star 
and is based on synthetic photometry in the expected (at the 
time) SDSS u'g'r'i'z' filt ers and an improved spectral energy 
distribution for this star (Fukugita et alJ[l996T) . This is used 
to calibrate a primary network of 158 stars observed by the 
USNO 40-inch at Flagstaff in Arizona, chosen to span a range 
in c olor, airmass, and RA, and distributed over the Northern 
sky dSmith e t al. 2002). Unfortunately, these stars saturate the 
SDSS survey telescope; the calibrations are therefore indi- 
rectly transferred via a 20 inch Photometric Telescope (PT) 
at Apache Point (APO), which observes these primary stan- 
dards, as well as 1520 41.5x41.5 arcmin secondary patches 
of sky. These patches are wha t finally calibrate t he data from 
the 2.5m survey telescope (Tuck eret alj 120061) . Note that 
these are three different filter systems, and not just realiza- 
tions of one system. In addition to conflating the absolute and 
relative calibrations, this indirect transfer of the calibration 
makes achieving < 1% calibrations via this method a chal- 
lenge, although it does re turn relative calibrations accurate to 
~ 2% dlvezic et alj|200 4). Note that these errors have natural 
scales of 2.5°/12 (the width of a camcol) and 2.5° (the width 
of a stripe) perpendicular to the scan direction. 

Before continuing, it is worth emphasizing that 2% calibra- 
tions for a wide field optical survey was unprecedented until 
very recently. However, motivated by the promise of future 
wide field surveys and the challenge of 1% photometry, we 
realized that the next step must short-circuit the above multi- 
stage calibration pipeline. The calibration algorithm we pro- 
pose here relies on repeat observations to constrain the pho- 
tometric calibrations. Unfortunately, in the standard survey 
strategy, the only significant repeat observations occur at the 



poles of the survey (where the great circles of stripes con- 
verge), and on the celestial equator, which is re-imaged every 
Fall (see Fig. |2j. While the Fall equatorial stripe has suffi- 
cient repeat observ ations to make precise photometry possible 
dlvezic et al.ll2007l) . the calibration for the bulk of the survey 
region would only be constrained at the survey poles, clearly 
undesirable. The only other natural overlaps occur when the 
beginning and ends of runs overlap each other along strips. 
While this does connect the survey from one pole to the other, 
most of the overlaps occur on the same CCD column, and so 
have limited utility since these are degenerate with flat fields. 

To address both these inadequacies, two additional sets of 
data were taken. The first were short scans that cross the nor- 
mal scan directions. Such oblique scans exist for most observ- 
ing years (Fall through Spring), and were taken to check for 
temporal variations of the flat fields. These are invaluable for 
constraining flat fields, since they compare each CCD column 
with every other. The other data were a grid of long scans, 
dubbed the "Apache Wheel", designed to connect every part 
of the survey with every other. Observing such a grid with the 
usual SDSS scanning speed would have required a significant 
expenditure of telescope time, adversely affecting the science 
goals of the survey. The compromise was to observe these 
data, at 7 times the normal scanning speed (i.e. with an ef- 
fective exposure time of ~ 8 seconds), and binning data into 
4x4 native camera pixels. Reducing these data r equired mod- 
ificati ons to the survey data reduction pipelines ( Lupt on et al.1 
l200l . and was done at Princeton (along with a re-reduction 
of regular survey data) as part of this calibration effort. The 
survey region we consider in this paper is in Fig. |2] with the 
greyscale encoding the number of repeat observations of dif- 
ferent regions of the sky. 

3. THE ALGORITHM 

3.1. The Photometric Model 

An introduction to photometric ca librations and photomet- 
ric standard systems may be found in Bessell (2005); we focus 
on the details relevant to this work below. Assuming linear- 
ity, the flux / of an object at Earth (above the atmosphere) is 
related to the detected flux Jadv 11 by 



/ = ICfADu ; 



(1) 



the problem of photometric calibration is to determine K-. The 
above equation is deceptively compact; K. depends on the ex- 
posure time, detector efficiency, filter responses, the telescope 
optical system, the optical path through the atmosphere, the 
spectral energy distribution of the objects in question, and all 
the variables that these in turn are sensitive to. Furthermore, 
Eq. Q] makes no reference to the units of / and /C; the prob- 
lem of determining the correct units is that of absolute pho- 
tometric calibration; we restrict our discussion below to the 
problem of relative calibrations. 

Since all the above terms affect the flux multiplicatively, it 
is convenient to work in log-space; the above effects become 
additive corrections. Converting fluxes to magnitudes (m = 
— 2.51og 10 /), Eq. [Ubecomes 

m = m ADU - 2.5 log 10 (/C) . (2) 

Expanding JC in terms of its various dependencies, we obtain 

-2.5\og w ()C) = a(t)-k(t)x + f(i,j;t) + ... , (3) 

17 An ADU, or Analog-Digital Unit is the digitization of the analog detec- 
tor output 




FIG . 2. — The sky coverage of the SDSS data used in this paper, shown in an equal area resolution 7 HEALPIX/HEALCART (Gorski et al. 1999; Finkbeiner 
2004) projection. The x scale covers RA 0° to 360°, while the y axis runs from DEC 90° to —90°. The grey scale denotes the mean number of observations of 
a star in a particular pixel. Note that we saturate at 5 observations, although on the equatorial (white) stripe, there are pixels with a mean number of observations 
as high as 15. The bulk of the survey data is in the North Galactic Cap, the prominent structure in the center of the image. The Equatorial stripe, imaged every 
Fall, is the white horizontal stripe halfway in the image. The approximately equally spaced vertical runs are examples of the Apache Wheel data. 



TABLE 1 
Flat Field seasons 



SDSS Run 


MJD 


Date 


Comments 


1 


51075 


19-Sep-1998 


Beginning of Survey 


205 


51115 


28-Oct-1998 




725 


51251 


13-Mar-1999 




941 


51433 


12-Sep-1999 




1231 


51606 


()3-Mar-2000 




1659 


51790 


03-Sep-2000 


After i2 gain change 


1869 


51865 


17-Nov-2000 


Vacuum leak in Dec 2000 


2121 


51960 


20-Feb-2001 


After vacuum fixed 


2166 


51980 


12-Mar-2001 




2504 


52144 


23-Aug-2001 


After summer shutdown 


3311 


52516 


30-Aug-2002 


After summer shutdown 


4069 


52872 


20-Aug-2003 


After summer shutdown 


4792 


53243 


26-Aug-2004 


After summer shutdown 


5528 


53609 


26-Aug-2005 


After summer shutdown 



NOTE. — The starting dates, and the corresponding first SDSS run for the 
fiat field seasons. 

where all terms are a function of time. The optical response 
of the telescope and detectors is the "a-term" a(t), while the 
detector fiat fields (in magnitudes) are f(i,j;t) where i, j rep- 
resent CCD coordinates. The atmospheric extinction is the 
product of the "k-term" k(t) and the airmass of the observa- 
tion, x. Note that this is a crude phenomenological model (it 
heuristically resembles a first order Taylor expansion), but is 
completely adequate for our purposes. We therefore defer a 
discussion of its limitations and potential extensions to Sec. [7] 
We now specialize to the SDSS; we calibrate each of the 
five filters individually, and assume that each of the six camera 
columns are independent, yielding an a-term and fiat field to 
be determined per CCD. We implicitly assume that the filter 
response for each of the six CCDs is identical (we return to 
this in Sec. 13. The k-term is however common to all camera 
columns and depends only on the filter. Also, since the SDSS 
observes by drift-scanning the sky, the fiat fields are no longer 
two-dimensional, but only depend on the CCD column and 
are represented by a 2048 element vector. This is complicated 
by the fact that some of the SDSS CCDs have two amplifiers, 
resulting in a discontinuity at the center of the flat field. We 



model this by assuming the fiat fields have the form, 

/(^) = /o(j)+0(j-lO24)A/ (4) 

where 6{x) is the Heaviside function, and A/ (hereafter, the 
"amp-jump") is the relative gain of the two amplifiers. Note 
that as written, fo is a continuous function of CCD column. 
Finally, we need to specify the time dependence of these 
quantities. The a-terms and amp-jumps are assumed to be 
constant during a night, and we simply specify these as piece- 
wise constant functions. 

It was also realized early in the survey (about 2001) that the 
flat fields were time-dependent, and appeared to be changing 
discontinuously over the summers when the camera was dis- 
assembled for maintenance. These changes are most likely as- 
sociated with changes in the surface chemistry of the CCDs. 
We therefore model the flat fields as being constant in time 
over a "flat field season", roughly the period between any 
maintenance of the camera. The boundaries, in MJD and 
SDSS run number, of these "seasons" are listed in Table Q] 
Ideally, one might have chosen an even finer time interval to 
test the constancy of the flat fields; however, the SDSS lacks 
sufficient oblique scan data to improve the time resolution. 
We note here that the standard practice of measuring flat fields 
from sky data does not work for the SDSS, due to scattered 
light in the camera. 

The time dependence of the k-terms at APO is more com- 
plicated, as the atmosphere (on average) gets more transparent 
as the night progresses, at the rate of ~ 1 mmag/hour (milli- 
magnitudes/hour) per unit airmass. We therefore model k(t) 
over the course of a night as 

dk 

k(t) = k + —(t-t ref ), (5) 

where t re f is a reference time 18 . Note that t in the above 
equation only runs over the course of a single night; k and 
dk/dt can (in principle) vary from night to night, and there is 

18 We adopt 0700 UT as t re r, corresponding to midnight Mountain Stan- 
dard Time. 
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TABLE 2 
Calibration Parameters 



Parameter 



Number 



Fit 



Comments 



a-terms 6 X 5 X n n i g b,t Yes 

k 5 X rinight Yes k-term at t = t r 

dk/dt 5 No 

Flatfields 6 X 5 X n sea son Yes 2048 element vector 

(iterative) 

Ampjumps 6 X 5 X n rurL No 



NOTE. — The parameters that make up the photometric model. The num- 
ber of parameters are functions of n n i g b.t (the number of nights), 

^■season 

(the number of flatfield seasons), n run (the number of rims), and the number 
of filters (5), and camera columns (6). Also listed is whether the parameter is 
fit for or not. 

no requirement on the continuity of k(t) across nights. Table[2] 
summarizes the parameters in our photometric model, whose 
final form is 



m = rriADU +a a ■ 



(6) 

with a, (3, and 7 indexing the appropriate a-term, k-term (and 
t re f), and flat field for the star in question. 

3.2. Solution 

Having specified the parameters of the photometric model, 
we now turn to the problem of determining them. It is natural 
to consider repeat observations of stars to constrain these pa- 
rameters 19 . Let us therefore consider n b s observations with 
observed instrumental magnitudes rriADU.j, of n s t a r unique 
stars with unknown true magnitudes mj. Note that n a b s is the 
number of observations of all stars, i.e. n Q b s = 2~Z"=i ar n « 
where n, is the number of times star i is observed. Using 
Eq. [6] we construct a x 2 likelihood function for the unknown 
magnitudes and photometric parameters, 



kp,(dk/dt)p,f^\ = E x 



(7) 



with 



x?= E 

jeO(i) 

(8) 

where j runs over the multiple observations, 0(i), of the i th 
star, a is the error in mj^ADU, an d k(t) is given by Eq.|5] We 
also assume that errors in observations are independent; this 
is not strictly true as atmospheric fluctuations temporally cor- 
relate different observations. One can generalize the above to 
take these correlations into account and, as we show below, 
our results are not biased by this assumption. Note that Eq. [7j 
has n b s known quantities, and n star + # (parameters) un- 
knowns. In general, the number of photometric parameters 
<C n s t a r, and n bs > 2n s t a r> implying that this is an overde- 
termined system. 

To proceed, we start by minimizing Eq. UJ with respect to 
mc, this yields, 



dm,: L3£ ° (l) 



mi—mj t ADU 



-q a( j ) +fc w) (t)x-/. T (j ) 
^ 



(9) 



While in principle, one could also consider galaxies, we restrict our 
discussion to stars to avoid subtleties of extended object photometry. 



which is trivially solved for m, to give, 



nv. 



E 



nij,ADU + a a {j) ~ kp(j)(t)x + 




•(10) 



As substituting the above result into Eq.|9]to solve for the cal- 
ibration parameters is algebraically unwieldy, we reorganize 
these results by making the following notational change. We 



arrange the unknown photometric parameters into an n pa 
ement vector p, 



el- 



kp 
(dk/dt)p 

fi 



(11) 



Then substituting Eq. [TOlinto Eq. [8] yields a matrix equation 
for x 2 , 



X 



= (Ap-ty'CT^Ap-b) 



' I hi r 



(12) 



matrix , and b is an n Q bs element 



where A is an n b s x 

vector, and v* represents the transpose of v. The errors are 
in the covariance matrix C which, in Eq. [8] is assumed to be 
diagonal (but can be generalized to include correlations be- 
tween different observations). For clarity, we explicitly write 
out the form of Ap — b for the case of a single star observed 
twice at airmass x\ and X2, and with errors o\ and 02, where 
only the a- and k-terms are unknown, 



1 
1 









~%2 



h 
h 



h 
h 



-X\h 

-x\h 



-X2I2 
-X2I2 



( at \ 

a.2 

h 
V J 

m^ADU - rni.ADuh - m 2 ,ADuh 
m 2 ,ADU - mi^ADuh - m 2 ,ADuh 



(13) 



where Ii is the normalized inverse variance, l. L — 
a^ 2 j Y^,j a J 2 - Each row of Ap — b has a simple interpre- 
tation as the difference between the magnitude of a particular 
observation of a star, and the inverse variance weighted mean 
magnitude of all observations of that star. Also, although A 
is a large matrix (~ 50, 000, 000 x 2000 for the SDSS), it is 
extremely sparse, and amenable to sparse matrix techniques. 

Obtaining the best fit photometric parameters simply in- 
volves minimizing Eq.Q~2] Although there are sev eral choices 
here, we proceed via the normal equations (e.g. JPress et alj 
fT992h . 

= A*C _1 Ap - A*C _1 b = . (14) 

dp 



The inverse curvature matrix, 



d 2 x 2 
dp l dp j 



= (A'C^A) 



(15) 



provides an estimate of the uncertainty in the recovered pa- 
rameters. Note that it is however not the covariance matrix of 
the parameters, since Eq. [12] was derived marginalizing over 
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the unknown magnitudes of all the stars. Furthermore, since 
the measurement errors do not account for temporal variations 
in the atmosphere, the "error" estimates from the curvature 
matrix may be significantly underestimated. 

We conclude by noting the similarities between the above 
and algorithms used for making maps of the CMB (e.g. 
lTegmark|[l997l) 20 . The SDSS runs are analogous to CMB 
scan patterns, while the magnitudes are equivalent to the tem- 
perature measurements. However, unlike the CMB, our prin- 
cipal goal is the calibration parameters, with the magnitudes 
of the stars being a secondary product 21 . 

3.3. Degeneracies and Priors 

Our choice of photometric parameters is non-minimal in 
that there exist degeneracies between them. These degenera- 
cies are of more than academic interest, as they make the nor- 
mal equations singular, and solutions of them unstable. We 
now discuss the source of these degeneracies, and how the 
resulting numerical instabilities can be tamed by the use of 
priors. 

• Zero point : As the above algorithm is based solely on 
magnitude differences, any overall additive shift of all 
the a-terms does not change \ 2 ■ Note that this is simply 
the problem of absolute calibration rephrased. 

• Disconnected Regions : This is a generalization of the 
previous case; the zero points of each disconnected re- 
gion of the survey can be individually changed, with- 
out changing \ 2 . Note that disconnected in this context 
refers regions with neither spatial nor temporal overlap 
(as we assume the photometric parameters are stable 
over the course of a night) with other parts of the sur- 
vey. 

• Zero point of flats : In Eq. [6] the zero point of the flat 
fields is degenerate with the a-terms; this degeneracy 
is trivially lifted by forcing the flat fields to have zero 
mean. 

• Constant Airmass : The photometric equation schemat- 
ically is ~ a — kx; therefore for data with little or no 
airmass variation, there is a degeneracy direction that 
keeps a — kx constant, while changing both the a- and 
k-terms. While this does not affect the calibration in 
regions where a — kx is constrained, extrapolating the 
a- and k-terms to regions with different airmasses can 
result in incorrect calibrations. 

There is a useful generalization of the above discussion; 
the inverse eigenvalues of the curvature matrix (Eq. [14] and 
[T5l ) are a measure of the error on the determination of lin- 
ear combinations of the photometric parameters (encoded by 
the corresponding eigenvectors). The degeneracies discussed 
above are characterized by eigenvalues ~ 0, which make the 
normal equations unstable. However, any badly constrained 
combinations (even if they are formally well determined) can 
amplify noise and un-modeled systematics in the data, poten- 
tially introducing errors when the calibrations are applied. We 
therefore identify all eigenvectors of photometric parameters 

20 This is not accidental, as this work was inspired by the techniques 
learned in CMB mapmaking. 

21 This results in the unusual situation of having ~ a few million nuisance 
parameters that must be marginalized over to obtain ~ a thousand parameters 
of interest. 



that are poorly constrained (i.e. those that could result in po- 
tential errors of > 1% ), and project these out; this renders 
the normal equations stable and they can be directly solved. 
Note that this introduces a tunable parameter to the solution - 
the eigenvalue threshold below which we project out modes. 
This threshold is chosen such that the final calibrations are 
insensitive to its exact value. 

Although projecting out poorly determined eigenvectors 
yields a minimal set of parameters well constrained by the 
data, we must add back in these "null" eigenvectors to get a 
solution in our original (and preferred) parameter space. We 
achieve this by introducing priors on the photometric param- 
eters, and then adjusting the values of all the photometric pa- 
rameters along the null vectors to best satisfy these priors. As- 
suming equally weighted Gaussian priors on the parameters p 
with a mean value po, this can be phrased as an auxiliary \ 2 
minimization, 

Xprior = IP' + V n „«<5x - p | 2 , (16) 

where p' is the solution of the normal equations, Eq. [14] 
and we have gathered the ridegen null eigenvectors into a 
n par x n degen matrix V nua . Varying 5x to minimize Xprior> 
we obtain our final solution for the photometric parameters, 

p = p' + V nuU 6x . (17) 

3.4. Implementation Details 

The above discussion described our calibration algorithm in 
generic terms, with minimal reference to survey specifics. We 
now discuss the details and approximations specific to imple- 
menting this algorithm for the SDSS. 

The first approximation involves determining the flat field 
vectors. As described, the flat field vectors are determined si- 
multaneously with the other photometric parameters. Doing 
so would however have approximately doubled the number 
of photometric parameters, and significantly complicated the 
degeneracies between the various parameters. We therefore 
chose an iterative scheme where the flat fields are held con- 
stant while the other parameters are determined. We then use 
the best fit solution to measure the magnitude differences be- 
tween multiple observations as a function of CCD column, 
and fit a flat field vector to these via a quadratic B-spline with 
17 uniformly spaced knots. As we show in the next section, 
this scheme rapidly converges to the true solution. In addi- 
tion, the SDSS photometric pipeline estimates the amp-jumps 
by requiring that the background be continuous across the am- 
plifiers. Instead of fitting to the amp-jumps, we simply hold 
them fixed to these values. 

The second approximation involves the k-terms and their 
time derivatives. The typical airmass variations over the 
course of a single night tend to be small, making the deter- 
mination of k-terms very degenerate with the a-terms, as dis- 
cussed in the previous section. The situation is even more 
degenerate for the time derivative of the k-terms. We fix these 
degeneracies by using priors for the k-terms, and fixing their 
time deriv atives to values es timated by the SDSS photometric 
telescope (Hogg et al. 2001). Table [5] summarizes which pa- 
rameters are fit in our implementation, while Table [3] lists the 
mean values for k and dk/dt. 

We must also specify the actual objects used for calibrat- 
ing. We restrict ourselves to objects that the SDSS classifies 
as stars, and use aperture (7.43 arcsecond radius) photometry 
to determine their magnitudes. The first choice sidesteps the 
subtleties of galaxy photometry, while aperture photometry 



Improved Photometric Calibration of SDSS Imaging Data 



7 



TABLE 3 



Filter 


Mag. 


f^star 




ko 


dk/dt 


a(dk/dt) 




Limit 


(xlO 6 ) 


(xlO 6 ) 




(xlO 2 ) 


(XlO 2 ) 


u 


18.5 


4.7 


14.6 


0.49 


- 1.2 


2.5 


9 


18.5 


9.3 


29.1 


0.17 


-0.7 


1.7 


r 


18.0 


11.7 


36.5 


0.10 


-1.0 


1.7 


i 


17.5 


11.5 


35.9 


0.06 


-1.2 


1.5 


z 


17.0 


11.6 


36.4 


0.06 


-2.2 


1.7 



NOTE. — The magnitude limits used to select stars for calibration for the 
5 SDSS filters, with the resulting number of unique stars, n B tar, and the total 
number of observations, n i a (in millions of stars). Also listed are the mean 
k-term ko (used as a prior), the mean time variation of the k-term, dk / dt (in 
mags/airmass/10 hours), and its scatter about the mean. The latter is used in 
our simulations to determine the step size for the random walk approximation 
to the atmospheric extinction. Note that we do not fit for the time variation of 
the k-term but simply use the values for the entire survey. 

avoids aliasing errors from the point spread function (PSF) es- 
timation into the calibration. The magnitude limits we use are 
in Table [3] along with the number of unique stars, and obser- 
vations. We choose not to make any color cuts on the stars to 
eliminate variable stars and quasars. These only add noise to 
any calibrations, but cannot bias the results; we therefore just 
use outlier rejection (3cr clipping) and iterate our algorithm to 
minimize such contamination. A significant advantage of this 
approach is that the calibration of the 5 SDSS filters is inde- 
pendent, allowing us to use colors of sub-populations of stars 
as external tests of the calibrations; this is discussed in detail 
in Sec. El 

Our algorithm does assume that the input data were taken 
under photometric conditions. We therefore, at the outset, 
eliminate all data taken under manifestly non-photometric 
conditions. As we discuss below, the algorithm does provide 
diagnostics of the photometricity of the data; we therefore it- 
erate the algorithm removing any remaining non-photometric 
data. 

Finally, there remains the issue of the absolute calibration 
of these data, or determining the five zero points for each fil- 
ter. Improving the absolute calibration is beyond the scope of 
this paper; we therefore determine the zero points by match- 
ing magnitudes on average to those obtained by the standard 
SDSS calibratio n pipeline. These are therefore essentially on 
an AB system dAbazaiian et al.ll2004. tied to the SP SS fun- 
damental standard, BD+17°4708 (lOke & Gunnll 19831) . 

4. SIMULATIONS 

Simulations serve the dual purpose of verifying the above 
algorithm and our implementation of it, as well as quantifying 
the level of residual systematics. We construct the simulations 
as follows : 

• We start with the actual catalog of stars observed by the 
SDSS, with the magnitude cuts described above. This 
ensures (by construction) that the pattern of overlaps in 
the simulations matches the observed data, essential to 
obtaining realistic results. 

• We simulate "true" magnitudes for each of the stars, 
using a power law distribution, where the normalization 
and slope are matched to their observed values. 

• Given an observation of the star, we then transform the 
magnitude into an observed instrumental magnitude, 
assuming values for the a- and k-terms and flat fields. 



TABLE 4 
Calibration Errors 



Filter 


(Am) 


cr 


o"3 


%(3<r) 


CO 


u 


-1.67 


13.38 


12.53 


0.85 


7.25 


9 


0.82 


7.79 


7.31 


0.72 


1.77 


r 


0.93 


7.81 


7.26 


0.81 


1.69 


i 


0.92 


6.84 


6.38 


0.75 


1.32 


z 


0.97 


8.06 


7.61 


0.68 


2.70 



NOTE. — A summary of the calibration errors for the five SDSS filters, 
as determined by simulations; all values are in mmag. (Am) is the mean 
of the difference between the estimated and true calibration value for each 
SDSS field, while a is the corresponding standard deviation, with 0-3 the 3<r 
clipped value, and %(3cr) the fraction (in percent) of 3<r outliers. Finally, 
(jo is the calibration error just from measurement noise (i.e. for a simulation 
with no unmodeled random component to the atmosphere). 

We simulate the time variation of the k-term by describ- 
ing k(t) — fc by a Gaussian random walk in time, with a 
drift in time given by dk / dt. The size of the steps is set 
by the observed value of the scatter in dk/dt (Table[3]l. 
Note that this random component attempts to model the 
correlations in time induced by the atmosphere, albeit 
by making the simplification that the spectrum of fluc- 
tuations is described by a Gaussian random walk. Non- 
photometric data is simulated by exactly the same pro- 
cess, although we arbitrarily increase the scatter in the 
random walk. 

• We add noise to the instrumental magnitudes by con- 
sidering the Poisson noise from both the object and the 
sky. Note that the Poisson fluctuations from the sky 
dominate the error budget for most of the objects. 

These simulated catalogs are structurally identical to the ac- 
tual data. We can therefore analyze them in exactly the same 
manner, and compare the derived parameters with those in- 
put, providing us with an end-to-end test of our pipeline. Fur- 
thermore, these simulations have exactly the same footprints, 
timestamps and overlap patterns as the real data, allowing us 
to estimate our final errors and explore parameter degenera- 
cies. 

4.1. Results 

Figs.|3]and|4]show the differences between the true and es- 
timated a- and k-terms, and flat field vectors for the r band 
in one of our simulations, analyzed identically to the real 
SDSS data. The flat field vectors are recovered with an error 

< 0.5%. The SDSS pipeline stores flat fields as scaled inte- 
ger arrays; the roundoff error from this is about an order of 
magnitude lower. The a- and k-terms are similarly correctly 
estimated on average, although there are significant misesti- 
mates for both. However, a striking feature of Fig. [3] is the 
similarity in the residuals for the a- and k-terms, reminiscent 
of the discussion of the degeneracies between the a- and k- 
terms in Sec. 13.31 This suggests comparing the estimated and 
true values of a — k(x) on a per- field basis, where (x) is the 
average airmass over a given field and filter; it is this combi- 
nation that determines the photometric calibration of a field. 

The results of this comparison are in Table |4] We start by 
noting that the calibrations are determined correctly (on av- 
erage) to ~ 0.1% or better, verifying both the algorithm and 
our implementation of it. The errors in the calibrations are 

< 1%, or 10 mmags for all the filters (except u, where they 
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FIG. 3. — The difference between the estimated and true a- and k-term for 
the r band in one of our simulations. There are approximately 6 a-terms that 
correspond to a given k-term, and the scales on the x axis are adjusted so that 
corresponding terms are aligned. Note that the estimated a- and k-terms are 
highly covariant. 



are slightly higher), suggesting that the SDSS can break the 
"sound barrier" of delivering 1% relative calibrations over the 
entire survey region. Catastrophic failures in the calibrations 
are also negligible, evidenced both from the near equality be- 
tween the sigma-clipped and total variances, and the almost 
Gaussian fraction of 3a outliers. Finally, we note that the er- 
rors in the calibrations are dominated by the unmodeled ran- 
dom fluctuations in the k-terms. Simulations with no random 
fluctuations achieve calibration errors of ~ 0.1%, suggest- 
ing that the SDSS calibration errors are therefore completely 
dominated by unmodeled behaviour in the k-terms. The ex- 
ception again is the u band where measurement noise is only 
a factor of ~ 2 smaller than the random noise in the atmo- 
sphere. 

The spatial distribution of the calibration errors is in Fig. [5] 
The calibrations are uniform across the whole survey area at 
the ~ 1% level, and are noticeably better at the survey poles 
where the number of overlap regions increases (see Fig. 0. 
Importantly, although there is spatial structure over individual 
SDSS runs (which is inevitable, given that we calibrate entire 
runs as atomic units), there are no coherent structures over the 
entire survey region. 

The above discussion assumes calibrations making the de- 
fault choices described in Sec. 13.41 We can use our sim- 
ulations to discuss the robustness of the algorithm to these 
choices below. For simplicity, we only consider the r band 
for these tests. 

• Magnitude Limits : As discussed above, the errors in 
the calibration are dominated by unmodeled systemat- 
ics in the atmosphere, and not measurement noise. We 
therefore expect the algorithm to be relatively insensi- 
tive to the choice of the magnitude limit. We explicitly 
verify this by re-calibrating after decreasing the magni- 
tude limit by 0.5 mag. Although this reduces the num- 
ber of stars and observations by 30%, the calibration 
errors are unaffected, as expected. 



FIG. 4. — The difference between the estimated and true flat field vectors for 
the r band of one of our simulations. Each line corresponds to a different flat 
field season. Since the mean of the flat fields are degenerate with the a-terms, 
we only plot the deviations about the mean. For clarity, only the flat field 
vectors for one camera column are plotted; the results for the other camera 
columns are similar. The errors in the flat field estimation are ~ 0.5% (peak 
to peak). 

• Apache Wheel data : As described in Sec.|2] the SDSS 
imaging data was supplemented by a grid of 4x4 binned 
data designed to improve the uniformity of the calibra- 
tion over the entire survey region. Calibrating the sur- 
vey without these data increases the calibration error to 
10.4 mmag (compared with the 7.8 mmag in Table |4)l, 
an increase of ~ 30%. Most of this increase is however 
driven by catastrophic failures; the 3er clipped variance 
only increases to 8.1 mmag, a more modest increase of 
10%. As expected, the Apache Wheel data better con- 
strain parts of the survey that were poorly connected, as 
they were designed to do. However, for regions already 
well constrained, the improvements are marginal. 

• dk/dt : Since we do not fit for a value of dk/dt, we 
must understand how errors in our assumed value of 
dk/dt propagate to the calibration. Fig. [6] shows the 
difference between calibrating a simulation assuming 
the correct value of dk/dt, and assuming dk/dt = 0. 
While the increase in the size of the calibration errors 
is small, the incorrect value of dk/dt introduces an 
overall tilt to the survey (in the figure, this is approx- 
imately 10 mmag). This tilt results from the fact that 
regions of similar RA are observed at approximately 
the same relative time in the night. The errors from an 
incorrect dk/dt therefore do not cancel, but accumu- 
late into a tilt, because we always observe the sky west 
to east. This is exacerbated by the fact that there is lit- 
tle data connecting the survey at the ends through the 
Galactic plane, and therefore no closed loops to pre- 
vent the appearance of such a tilt. This is the most se- 
rious systematic error in the calibration, and could af- 
fect any large scale statist ical measures. In fact, both 
iPadmanabhan et alJ (120071) and iBlake etall (120071) ob- 
serve excess clustering of photometrically selected lu- 
minous red galaxies at the very largest scales. We spec- 
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FIG. 5. — An image of the calibration errors for r band on the sky obtained for one of our simulations. The projection is the same as Fig. [2] but zoomed in on 
the Northern Galactic Cap of the SDSS. The greyscale saturates at magnitude errors of ±0.02 mag. 

compare our results, we describe both the internal consistency 
(Sec. 15. Il l and astrophysical tests (Sec. 15.41 ) we use to assess 
the photometric calibration. In addition, we als o ad dress the 
spatial structure of the calibration errors (Sec. 15. 2k as well 
as the photometric stability of the SDSS (Sec. [53b . Finally, 
we compare our calibrations with the currently public SDSS 
calibrations (Sec. 15.51 ). 

In what follows, we use "magnitude residual" to denote the 
difference between the (calibrated) magnitude of an observa- 
tion of a star and the mean magnitude of all observations of 
the star. 



5.1. Internal Consistency 

The first internal consistency test is the distribution of mag- 
nitude residuals. Since the scatter in the residuals also in- 
cludes measurement noise (a), it is more illuminating to con- 
sider x = (m — (m) )/o", if the measurement errors are a good 
estimate of the scatter in the residuals, \ should be Gaussian 
distributed with unit standard deviation. This is plotted for 
the stars used in the calibration, for the five filters, in Fig. [7] 
At the faint end, we observe that \ is distributed as expected, 
suggesting that the measurement noise is a good description 
of the scatter, and that calibration errors do not appreciably 
increase the scatter. The discrepancy at the bright end is due 
to a floor (a — O.Olmag added in quadrature) we impose on 
the magnitude residuals, to reflect the fact that the dominant 
error for these stars is no longer Poisson noise but possible 
systematics in the measurements. Note that calibration errors 




FIG. 6. — The difference in calibration between assuming dk/dt = 
and the true value. The tilt over the survey region is clearly apparent, and 
is approximately 10 mmag over the survey region. The greyscale goes from 
-0.01 mag to +0.005 mag. 

ulate that a tilt in the calibration could be a possible 
contaminant to the measurements on those scales. 

5. THE SDSS PHOTOMETRIC CALIBRATION 

Having described and verified our algorithm, we apply it to 
the SDSS imaging data. Since we do not have ground truth to 
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FIG . 7 . — The magnitude residuals weighted by their errors (x) as a function 
of apparent magnitude, for the five SDSS filters. The residual for a given 
observation is the difference between the observed magnitude and the mean 
magnitude averaged over all the observations of the star. The dotted lines 
show x = iL while the solid lines show the 16% and 84% contours; these 
should coincide with the \ = il lines if the scatter in the magnitudes is 
well described by the errors. The discrepancy at bright magnitudes is due to 
an enor floor we impose to down-weight the brightest stars. 

Flat #16: idFF-000725-r5.fit 
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FIG. 8. — An example of a flatfield vector from one r— band chip during 
Season 3. The grey scale and 25%, 50%, and 75% contours show the mag- 
nitude residuals as a function of CCD column, for all stars observed multiple 
times during that season. The smooth central (red) line shows the best fit 
(splined) flatfield vector. 



would only broaden the distribution of \- 

We also consider the magnitude residuals as a function of 
the CCD column, grouping the data by CCD and fiat field sea- 



FlG. 9. — Magnitude residuals as a function of time/field number for two 
example runs; all six camera columns are combined in these plots. The con- 
tours again show the 25%, 50%, and 75% levels. The red hatched regions 
mark periods of time independently known to be non-photometric from the 
SDSS photometricity monitors. Note that both these runs are on the multiply 
imaged Equatorial stripe, and therefore have lots of overlaps. However, for a 
large fraction of the data, the overlaps are considerably more sparse. 

sons; this is an estimate of the accuracy of our fiat field correc- 
tion. An example of these magnitude residuals as a function of 
the CCD column for camcol 5 in r— band is shown in Fig. [8] 
We don't correct for the flat field in this plot, to show the 
structure of the flat field itself. The r.m.s scatter in the magni- 
tude residuals about the derived vector is ~ 0.5% throughout 
the chip, although it increases at the edges of the CCD. Also, 
since we do not fit for amp-jumps but use the values derived 
from the photometric pipeline, the fiat fields adjust to correct 
for errors in the amp-jumps. Note that the errors in the amp- 
jump estimation are small (the true amp-jumps are usually a 
few tenths of a magnitude, while the errors are a few milli- 
magnitudes) that the splines have the necessary flexibility to 
adequately flatten the field. 

Finally, we plot the magnitude residuals, grouped by run, 
as a function of field number (and time); two examples are in 
Fig. [9] These plots are our primary diagnostic of the photo- 
metricity of the data. Photometric data have the mean residual 
scattered around zero, although often with coherent errors at 
the few millimagnitude level. By contrast, the residuals for 
unphotometric data show large excursions from zero, often at 
the ~ 10% level or greater. Most of these data have already 
been correctly flagged as being non-p hotometric by the SDSS 
photometricity monitors (Hogg et al. 2001) and have been ex- 
cluded from the solution. Any remaining non-photometric 
data is manually flagged as such, and removed in a second 
iteration of the calibration. For all the non-photometric data 
that overlaps photometric data, we can estimate an a-term per 
field that minimizes the residuals which determines the cal- 
ibration of those fields (these are still flagged as being non- 
photometric). 
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FIG. 10. — Examples of the spatial structure in the calibration errors for 
the r band, organized from left to right, and top to bottom in increasing order 
of their uncertainties. The top left mode is the best constrained, while the 
bottom right mode is the worst constrained. The middle row are examples 
of modes with typical errors. The modes are normalized such that the maxi- 
mum absolute error is 1 . Note that the worst constrained mode is the exactly 
degenerate overall zeropoint of the survey. The structures are similar for the 
other bands. 

5.2. Spatial Error Modes 

Since our goal in this paper is accurate relative calibration, 
it behooves us to understand the spatial structure in the cal- 
ibration errors. Our starting point is the curvature matrix, 
Eq. Q3] The eigenvectors of this matrix partition our basis 
of photometric parameters into uncorrelated linear combina- 
tions, whose uncertainties are given by the inverses of the cor- 
responding eigenvectors. An error in each photometric pa- 
rameter can be thought of as a pattern of errors on the sky, 
determined by the runs corresponding to that parameter. One 
can use this to project the eigenvectors (modes) of the curva- 
ture matrix on the sky. These then describe the spatial struc- 
ture of the calibration errors (Fig. ITOt . Note that projecting 
these modes on the sky destroys the linear independence of 
the modes; if desired, this can be restored by a straightfor- 
ward orthogonalization. 

The worst constrained mode is, as expected, the zero point 
of the calibration which is exactly degenerate. However, ex- 
amining the other poorly constrained modes (an example of 
which is the bottom left of Fig.flOb. we observe that there are 
no other such simple large scale modes, an indication of the 
fact the survey is well connected. At the other extreme are 
the best constrained modes. These are typically complicated 
combinations, and not surprisingly describe modes held to- 
gether by the grid of Apache Wheel data. More illuminating 
are examples of typical modes, two of which are in the middle 
row of Fig. [10] The most noticeable characteristic is the strip- 
ing along the scan direction. This simply reflects the fact that 
we calibrate camera columns individually, resulting in errors 
correlated in the scan direction. 

We do not fit for dk/dt, and so it does not ge t inc luded 
in the curvature matrix. However, we saw in Sec. 14. II that it 
resulted in a coherent tilt from one end of the survey to the 
other. 

5.3. Experimental Stability 

Our results for the photometric calibration can also be used 
to estimate the overall photometric stability of the SDSS cam- 
era, telescope and site. We estimate the camera stability by 
considering differences between a-terms as a function of time 
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FIG. 11. — [top] The difference between the a-terms of u band camera 
column 4 and 1; since the atmospheric corrections are common to both of 
these, this is a measure of the stability of the telescope + camera system. The 
drift in the camera during early data due to problems with the vacuum system 
(see Table [T} is clearly visible; the vertical dotted line at MJD 51960 marks 
when the vacuum system was fixed. Note that all of the changes are long 
term drifts; the system is stable on short (i.e. day) intervals as assumed in our 
model, [bottom] The drift in the relative (to camera column 1) zero points of 
camera columns 2 through 6, for all five filters in mmag/year, measured after 
MJD 51960. 

(for definiteness, we compute the a-terms relative to camera 
column 1); these differences are insensitive to any common 
mode effects (such as the atmosphere). An example is in 
Fig. QT| During the initial phases of the survey, we note that 
the camera was not very stable over long time periods, re- 
flecting various problems with the vacuum system flagged in 
Table [JJ However, over the past ~ 5 years, the camera has 
been extremely stable, as evidenced by an overall drift in the 
a-term differences of <^ 10 mmag/yr, for all the CCDs. 

One could also measure the combined stability, treating the 
camera, telescope, and site as a combined system. As the 
a- and k-terms are degenerate, we consider the combination 
a — k(x) every night, where we average the airmass over all 
the observations in a given night. This is plotted in Fig. [12] 
for the five SDSS filters. The most striking aspect of these 
data are the seasonal variations, seen as periodic oscillations 
in the data, at the ~ 10% level (except in the u band, where 
they are ~ 20%). Factoring out the seasonal variations, we 
find less than a 5% drift over the ~ 7 years considered here 
(again, except the u band, where the drift is ~ 10% over the 
same time period). 

We emphasize that all of the effects discussed here are long- 
term effects and do not affect the quality of the calibration 
which only assumes stability over a night for the a- and k- 
terms. 

5.4. Principal colors 

The above discussion has relied on a combination of simu- 
lations and internal consistency checks to assess the quality of 
the calibrations. While these provide essential perspectives, 
they have important disadvantages as well. Internal consis- 
tency checks are not independent of the calibration and might 
not flag deviations from the input model. Furthermore, these 
checks are local measurements, and do not provide informa- 
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FIG. 12. — The value of a — for camera column 1 as a function of 
MJD, and filter; (x) is the mean airmass of all the observations in a given 
night. This combination is insensitive to degeneracies between the a- and 
k-terms, and measures the overall photometric stability of the SDSS camera, 
telescope and site. The seasonal variations in these data are clearly apparent, 
as is the fact that the mirror is aluminized every summer. 

tion about large scale systematics problems. While simula- 
tions fill that gap, they are limited by the input model used. 
Astrophysical tests complement the above by providing large 
scale, independent verification, and are ultimately limited by 
astrophysical uncertainties. 

The majority (>~ 98%) of the stars detected by the SPSS 
are o n the main sequence (Finlat or" et alJ l2000l : iHelmi et alJ 
2003), and lie on one-dimensional manifolds (the "stellar lo- 
cus") on color-color diagrams. This suggests using the posi- 
tion of the stellar locus as a diagnostic of calibration errors 
dlvezic et al.ll2004t) . While there are a number of morpholog- 
ical fe atures one could use as a marker, we follow the discus- 
sion in llvezic et al.l ((2004) and use the "principal colors" that 
define directions perpe ndicular to the ste llar locus. We con- 
sider four such colors dlvezic et alj|2004l) : s (perpendicular 
to the blue part of the locus in u — g vs. g — r plane), w (the 
blue part in g — r vs. i — i), x (the red part in g — r vs. r — i), 
and y (the red part in r — i vs. i — z): 

s = -0.249u + 0.794 5 - 0.555r + 0.234 
w = -0.227 5 + 0.792r - 0.567i + 0.050 
x = 0.707.g- 0.707r- 0.988 
y = -0.270r + 0.800i - 0.534z + 0.054 . (18) 

We correct all magnitudes with the Schlege l et al.l d 1998b es- 
timates of extinction (except immediately below), but do not 
attempt any correction for stars not completely behind all the 
dust. Since we calibrate each band separately, and apply no 
color cuts to select the stars used, the above principal color 
diagnostics provide a completely independent verification of 
the calibration. 
Fig-QSplots these on a /i-stripe projection; the ^-direction 



is the coordinate along the scan direction /i dPier et alJ l2003), 
while the y coordinate is given by 

y = 12 (stripe) + 2(camcol) — 2 , strip = S 
= 12 (stripe) + 2(camcol) - 1 , strip = N , (19) 

where stripe, strip and camcol is the SDSS stripe number, 
whether it is a northern or southern strip, and the camera 
column respectively. This lays out each camera column as 
a row, respecting the interleaved structure of the strips within 
a stripe. The advantage of this projection is that calibration 
errors appear principally as stripes in the /i direction, while 
Galactic structure appears as irregular structures localized in 
/i, making it easier to separate the two. For the purpose of this 
plot, we use colors not corrected for extinction to highlight 
the Galactic structure. We simply exclude the small fraction 
of data not on survey stripes for the purposes of this analysis. 

We note that there is little visual evidence for any striping 
over the entire survey region in s, w, and x. Reddening from 
Galactic dust is clearly visible in s and x, which are colors 
nearly parallel to the reddening vector. The y map, on the 
other hand, does appear to show striping, with a periodicity 
on the SDSS stripe scale. In order to quantify this effect, we 
plot the average principal color per camera column (distin- 
guishing between northern and southern strips) in Fig. [T4l As 
anticipated from the 2D maps, the s, w, and x colors are uni- 
form at the 0.5% (peak to peak) level, whereas camera column 
2 is offset in y at 0.7% (peak to peak). It is unlikely that this 
is an artifact of the calibration process which treats all cam- 
era columns identically. Since y is the only color to use the 
z band, we speculate that this could be caused by the known 
variations in the z-band filter responses. However, as this ef- 
fect is of the same order as other systematics present in the 
calibration (and below our target of 1%), we simply caution 
the reader about this systematic in this paper, and defer its 
resolution to future work. 

The above has focused on large scale systematics; Fig. Q3] 
shows examples of the principal colors as a function of CCD 
column averaging over a random sample of runs in a flat field 
season. The deviations from a constant color are $1% for all 
colors, and < 0.5% for s and w, consistent with our estimates 
from simulations. We also observe errors in the amp-jump 
determination at the ~ 0.5% level, similar to Fig. [8] 

5.5. Comparison with previous results 

We conclude this section by comparing the calibrations pre- 
sented h ere with those publicly available as part of Data Re- 
lease 4 dAdelman-McCarthv et al.l [20061 DR4) 22 . Fig. [16] 
shows the difference between the aperture magnitudes of DR4 
and those derived in this paper, for all stars with r band mag- 
nitude less than 18. The magnitudes agree on average by con- 
struction, as the zero points were determined by matching to 
the public calibration. Furthermore, the scatter is approxi- 
mately 2% (r.m.s) for griz, and 3% (r.m.s) in u, consistent 
with the published uncertainties. The Data Release errors are 
therefore dominated by the Photometric Telescope (PT) based 
calibration method. 

Fig. [T7] plots these differences in the /i-stripe projection 
introduced previously. Since the standard SDSS calibration 
does not attempt to explicitly control relative calibration er- 
rors, the striping in the figure is not surprising. Note that the 
errors are correlated in the /i direction as expected, but also 

22 http : / /www . sdss . org/DR4 
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FIG. 13. — The spatial variation in the s (top left), w (top right), x (bottom left), and y (bottom right) principal stellar colors. The projection is a fi— stripe 
projection, with the x coordinate measuring /i (the coordinate along the scan direction). The y coordinate is the SDSS stripes, with each row as one of the 12 
camera columns that define the stripe. We have restricted ourselves to data on survey stripes between 9 and 44, corresponding to most of the North Galactic cap, 
in this plot. Note that the aspect ratio in the /i direction is significantly compressed. 



across camera columns. The latter arises from the fact that 
the calibration patches are 40' wide and span three camera 
columns, thereby correlating their calibrations. 

Finally, Fig.[l8]plots the differences in the DR4 flat fields, 
and those determined in this paper, for an example flat field 
season. The errors in the flat fields are both higher than 
the quoted uncertainties and appear to have long wavelength 
power. We speculate that these result from the method used 
to determine the flux response of the CCDs, which aliases 
flat field errors in the Photometric Telescope into the final flat 
fields. This aliasing is mitigated by using the average of g,r, 
and i, instead of any of those bands individually; this does 
not however eliminate the problem. Formally, these errors are 



~ 1% (r.m.s.), but are highly correlated, both spatially and in 
color. 

6. PUBLIC DATA RELEASE 

The calibrations (dubbed "ubercalibration") described in 
this work have been made public with the SDSS Data Release 
6 23 , and will be updated for the subsequent data releases. The 
SDSS Catalog Archive Server has recalibrated versions of the 
most popular magnitudes, as well as correction terms that can 
be applied to other magnitudes. We refer the reader to the 
documentation under "ubercalibration" on the SDSS data re- 

23 http : / /www . sdss . org/dr 6 
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FIG. 14. — The average principal colors measured for each of the 6 camera 
columns, for the north and south strips separately. As before, we restrict 
ourselves to stripes between 9 and 44. The mean color has been subtracted 
from each of the four curves; the means are -0.002, 0.004, 0.007 and 0.007 
for s, w, x and y respectively. The variations between camera columns are 
< 1% (peak to peak) for all colors. 
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FIG. 16. — Histograms of the differences in magnitudes between DR4 and 
this work, for stars with r-band magnitudes < 18.0. Shown are the distribu- 
tion of differences (normalized to have a maximum value of 1), as well as the 
cumulative distibution (dotted line). The vertical lines show the median of 
the distribution, which is < 0.001 mag for all five filters. 
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FIG. 15. — Examples of the principal colors as a function of CCD col- 
umn. The principal color and camera column used is noted in each panel. 
Each panel plots a color measured over a group of runs, chosen to be in the 
same Hatfield season (Table[T}. From 1 through 6, the runs used are (745, 
752,756),(1331,1345),(2190, 2299),(2566, 2662, 2883, 2886), (3560, 3830) 
and (4927, 5052) respectively. Deviations in the color from a constant (dotted 
line) indicate errors in the Hatfield determination. 




FIG . 17. — The difference in m agnitudes between DR4 and this work in the 
r band for the stars in Fig. 1161 plotted in the /^-stripe projection. The right 
panel zooms in to a region to highlight the structure in the calibration errors. 

lease websites for the most up-to-date information on these 
calibrations, and the available data formats. 

7. DISCUSSION 

We have presented a method for calibrating wide field 
imaging surveys using overlaps in observations, and applied 
it to the Sloan Digital Sky Survey imaging. Early versions 
of these calibrations have already been used for the creation 
of a number of auxiliary S PSS catalogs (e.g. Finkbein er et alj 
12004 iBlanton et al.ll200lh as well as a number of SPSS sci- 
entific publications (e.g.lTegmark et alJ2 004: Eise nstein et alJ 
120051: iPadmanabhan et alJl2007r.lTegmark et alj|2006l) . 

The principal features and results of this work are : 

• Relative vs. Absolute Calibrations: We explicitly sepa- 
rate the problem of relative calibrations from that of ab- 
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published SDSS flat fields are determined from the po- 
sition of the stellar locus. The use of stell ar flat fields 
mitigates all three of these (Manfr oTdl 1995b . Given suf- 
ficient observations and overlaps, one has sufficient S/N 
to map out the entire flat field with high precision. Fur- 
thermore, since one is using a realistic ensemble of stars 
by construction, biases due to differences between the 
flat field SED and the SED of a given object are mini- 
mized. Note that there are still potential biases for ob- 
jects of unusual color; these cannot however be treated 
in a general manner. 




FIG. 18. — The difference between the r band DR4 flat fields and those 
determined in this work, for SDSS runs between 4100 and 4400. These dif- 
ferences are traced by the difference in magnitudes as a function of CCD 
column, for the 6 camera columns. The solid line shows the median differ- 
ence, while the points are a subsampling of the individual measurements. We 
restrict ourselves to runs between 4100 and 4400 to select runs within a single 
flat field season. 

solute calibration. The problem of absolute calibrations 
then reduces to determining one zeropoint per filter for 
the entire survey (however, see the caveat on spectral 
energy distributions below), and does not alias into spa- 
tial variations in the calibration. This allows us to better 
control and quantify the errors in the relative calibra- 
tion. 

• Simulations : We emphasize the utility of simulations, 
both to validate pipelines, and to quantify the structure 
in the calibration errors. Simple analytical estimates 
are insufficient to characterize the errors, while astro- 
physical estimates are limited by their intrinsic scatter. 
Developing realistic simulations for the next generation 
of surveys must be an essential part of any calibration 
pipeline. Such simulations also are invaluable in de- 
termining the observing strategy {before any data are 
actually taken) that yields the desired calibration accu- 
racy. 

• Stellar Flat Fields : The problems of flat fielding wide 
field-of-view instruments, namely (i) non-uniform il- 
lumination for dome flats, (ii) spatial gradients and 
scattered light for sky flats and (iii) mismatched 
spectral energy distributions (SEDs) have been dis- 
cussed extensively in the literature (e.g. M anfroid 19951: 
Chromev & Hasselbacherlll996tlMagnier & Cuillandrd 
2004; Stu bbs & Tonry||2006l) . In particular, initial at- 
tempts to use sky flats for the SDSS resulted in errors 
of 5% in the r band and as bad as 20% in the u band, 
due to scattered light within the instrument. These were 
therefore never used for the public data; instead the 



• 1% Relative Calibrations/Spatial Error Modes : Our 
recalibrated SDSS imaging data has relative calibration 
errors, determined from simulations, of ~ 13, 8, 8, 7, 
and 8 mmag in ugriz respectively. We however do de- 
tect systematics not modeled in our simulations at the 
~ 0.5% level, suggesting a conservative estimate of 1% 
errors in griz and 2% in u. In addition, we are able to 
characterize the spatial structure of the errors as a com- 
bination of error modes. Most of these modes show lit- 
tle coherent spatial structure. The most significant spa- 
tial structure results from misestimating the time varia- 
tion of the k-terms, which introduces a tilt into the sur- 
vey. 

Throughout this paper, we ignored a number of sub- 
dominant systematic effects. We briefly discuss these below, 
both to document their existence as well as to alert future sur- 
veys of potential pitfalls. 

• Spectral Energy Distributions: When interpreting our 
magnitudes as absolute, our algorithm implicitly as- 
sumes that all objects have the same spectral energy 
distribution. The median r — i color of stars used in the 
calibration is ~ 0.2, and we expect the calibrations to 
be accurate (at the stated levels) for objects with colors 
not very different from these stars. This will however 
not be true for objects with unusual SEDs (e.g. SNe). 
In these cases, one must integrate the SED over the sys- 
tem response, in order to get an absolute, calibrated flux 
(inerg cm~ 2 s _1 ). 

• Filter Curves: We assume that the six copies of each 
filter are identical, and do not attempt any color correc- 
tions between the six camcols. We verified the validity 
of this assumption by generating synth etic ugriz pho- 
tomet ry for the Gunn-Stryker spectra dGunn & Strvkerl 
1983) for each of the six camcols, using the individu- 
ally measured filter curves. For stars with median r — i 
color close to ~ 0.2, the differences between the vari- 
ous camcols was better than 1% for griz and ~ 1% for 
the u — band. Of griz, the most drastic variation with 
color occurs for the z— band, with a 0.01 mag gradient 
between r — i = and 1 seen is almost all the camcols. 
Gradients of a similar magnitude are also seen for g2, 
g4 and r3. 

• Absolute Calibrations: We ignored the issues of de- 
termining the absolute calibration (i.e. the five zero- 
points) of the SDSS system, choosing instead to have it 
agree with the published SDSS magnitudes. In partic- 
ular, any corrections to put the SDSS system on to the 
AB system also apply in our case. 
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We conclude by discussing how to extend the program pre- 
sented here to the next generation of imaging surveys. Our 
starting point will be the second distinction made in Sec.Q]- 
separating the telescope and the atmosphere explicitly in the 
calibrations. It is relatively straightforward to adapt the al- 
gorithm presented here to use high precision measurements 
of the telescope response functions as a starting point; these 
would be analogous to the priors already considered here. 

Understanding atmospheric variations is an important step 
towards 1 % photometry; unmodeled variations are responsi- 
ble for almost all our calibration error budget. These trans- 
parenc y variations are domina ted by three well-studied pro- 
cesses (lHayes & Latham! 1975): Rayleigh scattering, molecu- 
lar absorption by ozone (dominant in the UV) and water va- 
por (dominant in the red and IR), and aerosol scattering. Of 
these, Rayleigh scattering is best understood, and is well de- 
termined by the local atmospheric pressure. While absorp- 
tion and aerosol scattering are well understood in an average 
sense, their time variation is significant. Tracking these would 
therefore require continuous monitoring of the atmosphere, 
plus detailed atmospheric models. The payback for doing so 
would be a dramatic reduction in calibration errors. 

The algorithm we propose in this paper demonstrates that 
1 % relative photometry is achievable by the current genera- 
tion of wide field imaging surveys. The challenge for the next 
generation of surveys is to break through the 1 % barrier. 
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